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Abstract 

The histogram reweighting technique, widely used to analyze Monte Carlo data, is shown to 
be applicable to dynamic properties obtained from Molecular Dynamics simulations. The theory 
presented here is based on the fact that the correlation functions in systems in thermodynamic 
equilibrium are averages over initial conditions of functions of the trajectory of the system in 
phase-space, the latter depending on the volume, the total number of particles and the classical 
Hamiltonian. Thus, the well-known histogram reweighting method can almost straightforwardly 
be applied to reconstruct the probability distribution of initial states at different thermodynamic 
conditions, without extra computational effort. Correlation functions and transport coefficients 
are obtained with this method from few simulation data sets. 
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The histogram reweigi 



;hting (HR) technique 



is, at present times, widely used in the 



determination of properties of thermodynamic equilibrium of many-body systems from a 
reduced amount of Monte Carlo data. It is recognized that the use of HR techniques increases 
the accuracy in the evaluation of thermodynamic properties and reduces the computer time 
needed for precise estimates of these properties. In this letter we address the problem of 
how the HR method can be extended to equilibrium dynamic properties, like correlation 
functions and transport coefficients, from molecular dynamics simulations. 

In a many-particle system in contact with a reservoir keeping constant a given set of 
intensive thermodynamic variables, the fluctuations in the conjugate extensive variables 
permit the system to explore thermodynamic conditions that differ from those fixed by the 
reservoir. The a priori knowledge of the functional form for the probability distribution of 
these fluctuations then permit to reconstruct, reweight, the histogram of the density of states 
of the system from a collection of Monte Carlo data sets. Then, the desired thermodynamic 
averages in different thermodynamic conditions can be evaluated. In the following, we will 
discuss how the same principle applies for dynamic properties of systems in thermodynamic 
equilibrium and show some results of the application of the method. 

The dynamics of a classical N-body physical system in phase-space is described by the 
Liouville equation 



where P(T,t) is the full phase-space probability distribution, H({pi}, {r*j}) is the Hamil- 
tonian of the system, while r*j and pi are, respectively, the position and momentum of the 
i t?l -particle. We are implicitly assuming that there are no external time- dependent force fields 
and that positions and momenta are generalized coordinates, for simplicity. We denote by T 
a given point in phase-space, i.e., to the ensemble of positions and momenta of all the parti- 
cles, {pi}, £(T) is the Liouville operator which governs the evolution of the probability 
distribution in phase-space. If we denote by A(T) and B(T) given phase-space functions, 
we can evaluate the correlation function from the two-event probability distribution, which 
can be recast in terms of the single-event probability and the conditional probability of two 




iCP(T,t) 



(1) 
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events, to yield 




(2) 



In this expression, P(T', 0) is the probability distribution of the system at the initial time, 
while the conditional probability P(T, t\V, 0) stands for the probability of finding the system 
in the state T at the time t, given that it was in the state r' at t — 0. If the system is in ther- 
modynamic equilibrium then P(T',0) = P(T',t) = P eq (T'). The form of the latter is known 
a priori and is characterized by the set of intensive variables kept fixed by the reservoir, 
which determine the actual thermodynamic state. For simplicity, we will only deal with the 
canonical ensemble, thus P eq (T') = exp(—H[T']/kT)/Z, with Z = J dT' exp(—H[T']/kT), k 
being the Boltzmann's constant. In eq. (0) P eq (T') is in fact the distribution of initial states 
of the subsequent dynamics of the system. It is crucial in our analysis to realize that P eq (T') 
contains all the required thermodynamic information of the system in the calculation of cor- 
relation functions. Effectively, the conditional probability P(T, t\V, 0) = exp(— itC)5(T — V) 
depends only on phase space variables. This conditional probability can in addition be writ- 
ten in terms of a trajectory in phase-space T[t; V], i.e. P(T, t\T', 0) = 5(T — T[t; T']), due to 
the fact that the Liouville operator simply introduces a displacement in the T coordinates, 
according to the Hamilton's equations of motion for a deterministic conservative dynamics. 
In these expressions, r[t; T'} is the set of positions and momenta of all the particles at the 
time t, which depend on the actual time as well as on the initial positions and momenta, T'. 
Hence, after integration over F the correlation function takes the form 



Eq. (JHJ) shows that (A(t)B(0)) depends on an equilibrium average of the function of the 
trajectory ^[tjT'] = A(T[t, T'])B(T') over the initial conditions. The thermodynamic in- 
formation is therefore contained only in the specification of these initial conditions. The 
function \I/ [t; T'] progresses along a trajectory according to the action of the Hamiltonian 
H[T]. If the Hamiltonian is time-independent, then the energy is therefore a constant of 
motion during the trajectory and so they are V and N. As a consequence, the evaluation 
of correlation functions corresponds to an ensemble average of trajectories characterized by 
constant values of N, V, and E, weighted with the probability distribution of the initial 
states corresponding to the actual thermodynamic ensemble. 




(3) 
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Therefore, from the simulation point of view, let us assume that we perform a set of 
different NVE molecular dynamics runs from specified initial conditions, randomly selected 
from a given thermodynamic ensemble. The values of the function ^[tjT'] are recorded, 
together with the values of the extensive parameters that characterize the statistical weight 
of the initial conditions T' of the given trajectory in the chosen ensemble. In the particular 
case of the canonical ensemble, the thermodynamic state is then characterized by N, V and 
T. Thus, we will record \l/[t; T'], together with E = H[T'] for every sample, due to the fact 
that P eq {T') oc exp(— H[T']/kT), in this case. In fact, an usual constant temperature MD 
simulation can be performed by using standard methods. The thermostat (a weak coupling 
with a bath) can be switched on, to appropriately set the initial state according to the suited 
canonical distribution, and off, during the measurement of the trajectory function 
Then, in a single simulation, the correlation function is obtained as usual from the simulated 
data, i.e. 

(A(t)B(0)) = J dV P eq (T'Mt;T'} ~ I",] (4) 

i=l 

In this expression, 1^ is the i th sample obtained in the simulation and A is the total number 
of samples. It is then clear that ^[tjT'] is a phase-space function that depends on one 
additional parameter, i.e. the time. Therefore, \& is suitable for the application of the HR 
technique provided that the sampling of the initial state distributions are propertly done. 

Effectively, in the standard application of the HR method, one would aim at reconstruct- 
ing the order parameter probability distribution 

P(*(t),T) = J dT'P eq (T',T)5^(t)-^[t;T'}) (5) 

at given thermodynamic conditions, from MD estimates of Since we are discussing the 
Canonical ensemble case, we have explicitly indicated the temperature T in eq. (J^J) and in 
what follows, for clarity reasons. Hence, from a set a — 1 ... N of independent molecular 
dynamics simulations, characterized each one by its temperature, T a , following the standard 
application of the HR method, one would write j^J 



Ea=idn a (t)P^(t),T ) 

Ea=l^aPa(^(t),T a 



p(*(t),T)d*(t) = ^ A „; T ,"; t\ (6) 



where A a is the total number of samples in the a th simulation and dn a (t) stands for the 
number of samples with a value of the variable in the range (^(t), *&(t) + d$fr(t)) in the a th 
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simulation. Moreover, P a (^>(t),T a ) is the order-parameter probability distribution in the 
a th data set. Finally, the desired correlation function is obtained from the reconstructed 
probability 

(A(t)B(0)) = [d$(t)*(t)P(*(t),T)~ 



Eq. ( j7j l. nevertheless, poses the difficulty of constructing the histograms for the variable 
ty(t). However, one can formally avoid the computation of histograms by inverting the 
order of summation to recover the sum over states|5j]. One then finally writes 

(A(t)B(0)) ( 3 T p^r T J, ) (8) 



a=l r' 



where P a (T') = exp(— H[T']/kT Q )/Z a , with Z a = J dT' exp(—H[T']/kT a ) and similarly with 
P(T',T), for the example discussed in this letter. As in the usual application of the HR 
method, the set of partition functions Z a has to be obtained following the standard proce- 
dure, which customary implies iteratively solving the corresponding integrals 0,0. Eq. (JEJ) 
is the most important result of this letter. On one hand, permits to obtain the value of a 
correlation function at a temperature T from a set of simulations made at temperatures T a . 
The expression is not limited to the canonical ensemble, but other ensembles can be analo- 
gously treated. On the other hand, the use of generalized Einstein relations or Green-Kubo 
expressions permit the evaluation of transport coefficients from the correlation functions. 

To illustrate the theory developed here, we present results for some correlation functions 
and transport coefficients of supercritical Argon near the critical point. The simulated 
system consists of iV = 256 Lennard- Jones particles with parameters e = 0.9961K J/ mol 
and a = 3.405A, in a cubic box of size L = 32.15646A so that the density of the system 
is p c = 510.7136% /m 3 , which corresponds to the critical isochore of the model used, being 



the critical temperature T c = 157.6577^ [6J. We have performed 10 MD simulations at 
temperatures ranging from T\ = 158.157-fT, slightly above the critical point, up to T2 = 
202.657i^. The simulations are performed thermostating the system via a Nose-Hoover 
procedure j^r] during hps., corresponding to 5000 time-steps. After this period, the energy of 
the initial state is recorded and a constant N, V, E integration is carried out for additional 
lbps., to obtain a sample of the suited function ^>[t; V] along the simulated trajectory. The 



process is repeated up to 2500 samples are obtained per simulation. Then the 10 MD data 
sets are used in eq. (JHJ) to interpolate the values of the desired correlation function at 50 
different evenly distributed temperatures in the temperature range between T\ and T^. In 
Fig^we can observe the velocity autocorrelation function for some selected temperatures 
in the analyzed domain. Upon integration of the velocity autocorrelation function, the 
corresponding Green-Kubo formula yields the diffusion coefficient of the system obtained 
as a continuous function of the temperature. The results are given in Fig. |21 
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FIG. 1: Velocity autocorrelation function (VACF) for 10 temperatures in the range 158.157-fT — 
202.657-ftT, selected every 5K, as obtained by means of the HR technique. The inner plot shows 
the initial behavior of the VACF. 



The evaluation of the viscosity of the system as a continuous function of the temperature 
has been carried out through an Einstein relation that involves the time-integral of the stress 
tensor 



V 



IV, d 
hm 



(9) 



20 k B T t-oc dt 

Here, a and (3 are indexes running over the three Cartesian coordinates, V is the volume, 
and AP a p(t) denotes the displacement of the elements of the pressure tensor P a/3 



AP a p (t) = ^d}- \P aP (r) + P Pa (r) - ^5 aP ^ Ppp (r)j 



(10) 



where 
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FIG. 2: Self-diffusion coefficient. The circles are the values of D obtained from the 10 MD 
data sets through the Green-Kubo formalism. The solid line are the estimates obtained from the 
application of the HR method using these data sets. Squares are independent MD simulations at 
T = T c , 169.6571T and 181.657iT, not included in the HR. The incept shows the behavior of D 
near the critical point. 



is the microscopic expression of the stress tensor |7j. In eq. (fTT|) . p a i is the a-component of 
the momentum of particle i, f a ij is the a-component of the force exerted on particle i by 
particle j, and r^j is the /^-component of the particle-particle vector, r^- = rj — iy Here, 
the trajectory function is AP a/3 (t), according to eq. (JTUJ). After the proper averages at a 
desired temperature are made using eq. (jHJ), then eq. (pj) yields the corresponding value of the 
viscosity coefficient at this temperature. Repeating this procedure for each temperature, one 
obtains a continuous curve, as can be seen in Fig. El The agreement between the simulated 
and empirical dataj^] is very good. In the strict vicinity of the critical point the results are 
obviously affected by finite-size effects that impeach the weak divergence of this coefficient 
at the critical point to be reproduced by our simple simulations For the same reason, 
no anomaly in the self-diffusion coefficient is observed [111]. A deep analysis of the dynamic 
behavior of the system at the critical point is beyond the scope of this letter. However, we 
want to stress that the method developed here might be very useful in the analysis of the 
dynamic critical behavior when applied in combination with finite-size scaling procedures. 
Finally, we want to point out that the HR method implemented here is sensitive to the 
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FIG. 3: Viscosity coefficient for Argon at supercritical conditions. The simbols have the same 
meaning as in Fig. EJ Dashed line corresponds to the experimental values 0. The inner plot shows 
the behavior of the shear viscosity close to the critical point. 

quality of the energy fluctuations produced by the thermostating method. In particular, we 
have observed that the energy fluctuations obtained from the method of velocity rescaling 
of Berendsen|l2| significantly differ from the Boltzmann distribution that we have assumed 
in the application of eq. (JHJ). As a consequence, the predictions near and out of the edges 
of the temperature range simulated, almost entirely constructed on these fluctuations, yield 
rather poor estimates of the physical quantities evaluated. 

Therefore, data on dynamic properties from MD simulations can be combined to extract 
dynamic information on the system at different thermodynamic conditions. We have shown 
that the method is suitable for the evaluation of correlation functions as well as the de- 
termination of transport coefficients through both Green-Kubo and Einstein relations. In 
particular, the agreement between the experimental values of the viscosity coefficient and 
the HR predictions from simulations is very good. In fact, as expected, the HR data shows a 
smoother behavior than the single simulation values of this coefficient since a larger amount 
of data is used in the interpolation than in the single simulation evaluation. The method 
is not limited to the Canonical ensemble but any other ensemble can be used if the proper 
probability distribution for the fluctuating extensive variable is a priori known. 
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